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We obtained a simple theoretical unified parameter for the characterization of rotational 
population propagation of diatomic molecules in a periodic train of resonant optical pulses. 
The parameter comprises the peak intensity and interval between the pulses, and the level 
energies of the initial and final rotational states of the molecule. Using the unified pa¬ 
rameter, we can predict the upper and lower boundaries of probability localization on the 
rotational level network, including the effect of centrifugal distortion. The unified parameter 
was tentatively derived from an analytical expression obtained by performing rotating-wave 
approximation and spectral decomposition of the time-dependent Schrodinger equation un¬ 
der an assumption of time-order invariance. The validity of the parameter was confirmed 
by comparison with numerical simulations for isotope-selective rotational excitation of KCl 
molecules. 
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I. INTRODUCTION 


The rotation of small molecules has been studied from the viewpoint of non-adiabatic molecular 


alignment [l-lsj] or molecular orientation [4l4lfl||. Non-adiabatic molecular alignment has been exper¬ 
imentally realized using non-resonant broadband optical pulses [ll, Recently, a train of pulses, 
whose repetition intervals are synchronized with the classical rotational period of the molecules, 
have been shown to be efficient for the rotational excitation of diatomic molecules Q, Q. At 
present, because of the development of terahertz-wave technol ogy ll 7 1. molecular orientation 
by rotationally resonant terahertz pulses has also been studied |l8l42ll|. The elemental technology 
employed to create a train of pulses in the terahertz region has already been reported 2^. An 


experiment involving molecular orientation using a train of rotationally resonant pulses will be 
performed in the near future. 

Using a train of optical pulses, we can impart the isotope selectivity in rotational excitation, 
because the classical rotational period depends on the mass of the molecules 


23h 27I|. The isotope 


selective population displacement o 


larly considered to be practical 


■ rotationally resonant terahertz pulses by the train is particu- 


28l-l3Cll|. Because the consecutive rotational population transfer by 


the train of pulses simultaneously starts from all the rotational states in the molecule, the principle 
of isotope selection still works well even if the initial rotational energies of the molecular ensemble 
are thermally distributed, specifically at high temperatures. However, there are still many prob¬ 
lems with the practical application of isotope separation of optical pulses in the terahertz region 
by the train. Further theoretical studies and development of elemental technologies are required. 

In real diatomic molecules, the consecutive rotational excitation by the pulse train has a certain 
upper boundary even if the interval of pulses is synchronized with the rotational period 3ll-l33l|. 


For a single pulse excitation, the limitation ordinarily comes from the spectral bandwidth of the 
pulse. For a multi-pulse scheme, two additional mechanisms induce the localization of rotational 
distribution; centrifugal distortion and mismatching of the pulse interval. Figure 1 shows an 
example of the time-evolution of the rotational probability distribution of diatomic molecules in a 
periodic train of pulses. For the isotope in a synchronized pulse train (Fig. UKa)), the wave packet 
of rotational probability is reflected at a certain limit, and repeatedly excited and de-excited inside 
the region. This phenomenon is localization by centrifugal distortion, called Bloch oscillation by 
FloB et al [^. For another isotope (Fig. HKb)), the population distribution was localized between 
J = 4 and J = 6. This is localization by interval mismatching, called Anderson localization by 


Flofi et al 


33|. 
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FIG. 1. Examples of numerical simulations of the evolution of rotational distribution for J = 5 of 

(a) and (b), by irradiating the resonant pulse train synchronized with rotational period of ^®K^^C1. 

The color scale is the same in both (a) and (b). 

We can intuitively understand the localization in the train of optical pulses by regarding the 
interaction as the resonant transition induced by the optical frequency comb. The train of optical 
pulses whose interval is Tp composes the optical frequency comb whose frequency peaks are evenly 
spaced by 1/Tp. Ignoring the centrifugal distortion, the one-photon pure rotational transition 
frequencies of the diatomic molecules are also evenly spaced. Isotope selectivity in rotational 
excitation originates from the agreement of the peak frequencies between the optical frequency comb 
and the molecular rotational transitions. Mismatching of Tp removes the comb’s peak frequencies 
from the resonance, and the consecutive resonant transitions are immediately prevented. The effect 
of centrifugal distortion can be ignored around the ground state; however, this effect becomes 
critical in the rotationally excited states. The centrifugal distortion shifts the transition frequency 
in proportion to (J + 1)^. Such a shift cannot be compensated by the optical frequency comb as 
long as only Fourier-transform-limited pulses are used. 

A quantitative evaluation of the region of localization of rotational distribution is important in 
the planning of experiments. This evaluation can be conducted by calculating the evolution of the 
time-dependent Schrodinger equation system; however, the dependence of the region of localization 
on the intensity of the pulses is odd and quite unclear. Therefore, to obtain the upper and lower 
boundaries of localization, it is necessary to perform numerical simulations for each case. 

In this study, we propose an evaluation method for the localization of rotational distribution 
without performing numerical simulations. We considered an evenly spaced Fourier-transform- 
limited pulse train in the terahertz region. We create a mathematical model using approximations 
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to the time-dependent Schrodinger equation system for diatomic molecules, and tentatively derive 
an analytical solution by partly neglecting precision. We extract a critical unified parameter from 
the analytical solution and validate it through comparison with the result of numerical simulations. 
In Sec. 2, we begin with an introduction to the targeted physical system and theoretically derive 
the unified parameter. In Sec. 3, the method used for the numerical simulations is provided. 
In Sec. 4, we show the results of the numerical simulations to be compared with the theoretical 
predictions. In Sec. 5, we discuss the results and meaning of the theoretical analysis. In Sec. 6, 
we provide the conclusions to this study. Atomic units h = me = e = 1, h = 2 tt are used unless 
explicitly stated. 


II. MATHEMATICAL MODELING 

A. Physical system 


We consider a one-dimensional network comprising a series of rotational states in diatomic 
molecules, which are connected by resonant optical transitions. In the interaction picture, the 
time-dependent Schrodinger equation of the diatomic molecules in a linearly polarized electric field 
is given by 

= -s{t)[nj-v exp{{Ej - Ej_i) ■ it} ■ Cj-i{t) 

+ /UJ -expKEj - Ej+i) ■ it} ■ Cj+i{t)], (1) 


where Cj{t) is the complex amplitude of the rotational states J at time t, e{t) is the electric field, 
fij is the transition dipole moment from J to J + 1, and Ej is the rotational energy of state J. In 
the diatomic molecules, the rotational energy Ej and the transition frequencies from J to J -|- 1 
are given by 


Ej — 2TrB]^J(^J -|- 1) — 2tt -|- 1)^, 

(2) 

= {Ej+i - 5j)/27r, 


= 2By[{J + 1) “ ‘iDy(j + 1)^, 

( 3 ) 


where By[ and are the molecular spectroscopic constants for a rotational series. Generally, the 
ratio of these constants is very small {D^/Bm — 10“® — 10“^). The classical rotational period of 
the diatomic molecules is defined by 


Tu 


1 

25m’ 


( 4 ) 
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ignoring centrifugal distortion. The molecular transition dipole moment is expressed as 

= fJ-o 


(J + l)2 -M2 


(5) 


(2J + l)(2J + 3)’ 

where fiQ is the permanent dipole moment, and M is the magnetic quantum number of the rotational 
series. Because the electric field is supposed to be polarized linearly, we need to consider only one 
Af in each series. In our mathematical treatment, we approximate fij as 


Mo 

y 


( 6 ) 


This approximation provides an almost exact solution for Af = 0, and can be used for small M in 
absolute values. Using the expressions given here, eq.([T]) can be rewritten as 


^ -y^[exp{27rj/j_i ■ it} ■ Cj-i{t) 
+ exp{-27ri/j • it} ■ Cj+i{t)]. 


We define the train of optical pulses using a multiple cosine electric field as 

N-l 

1 + 2 cos{27r • 2Bf{j + l)t} 


eit) = ^ 

Mo 


j=0 


(7) 


( 8 ) 


where 7 is the scaled intensity of the electric field, is the frequency constant that corresponds to 
Bm of the molecules, and N is the maximum J state to be connected by optical transitions. Under 
a rough estimation, the value of 7 corresponds to increments of J in a single pulse interaction. An 
example of the time-profile of the train of pulses is shown in Figi2j The interval between the pulses 



FIG. 2. A typical time-profile of a train of optical pulses, calculated using eq.® with N=100. 


and the maximum amplitude of a pulse are given by 

T 

2a’ 


^ma,x — 


f 

7-Sf 

Mo 


(9) 


(2A^ + 1). 


( 10 ) 
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B. Rotating-wave approximation 


From here, we perform a mathematical treatment to obtain an analytical expression for the 
rotational probability distribution. First, we apply the rotating-wave approximation, in which 


only the slowest oscillating terms are used and other quickly oscillating terms are ignored. We 


previously used this approximation to find the analogy between consecutive rotational excitation 


(CTQW) 


in a pulse train and continuous-time quantum walk 


. Under a fully resonant condition 


(Tp = Ty[-,D^ = 0 ), the evolution of the rotational distribution was observed to agree with the 
evolution of the CTQW. In our numerical calculations, under the conditions in which localization 
occurs, the rotating-wave approximation derives almost exact numerical results for a relatively 
weak electric field (7 < 10 ). 

Substituting eq.® into eq.©, and applying the rotating-wave approximation, we obtain the 
Schrodinger equation as 


= --^-^[exp{(i^j_i - dvrBfJ) • it} • Cj-x{t) 

-Fexp{(-i/j+i 47 rSf(J -h 1)) • it) ■ Cj+i{t)]. 


( 11 ) 


Without the approximation, eq. m involves many terms that oscillate with various frequencies 
in the exponential. The slowest oscillating terms remain in the exponential by considering the 
probable experimental conditions, that is, Bf = Bm and -C Bf. 


C. Decomposition of interaction matrix 


Using matrix representation, the interaction Hamiltonian of eq. m is given by. 


0 e-“o** 


0 


QOtQZt Q g OL^H 





( 12 ) 


Q Q—O-N-lit 


0 


gajv-lit g 


where aj is defined by 


aj = 27r{2AH(J + 1) - 4T»^(J -f 1)^}, 


(13) 
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where AB = Bm — Bf. The interaction matrix eq. (11211 is a time-dependent tridiagonal matrix, and 
is actually inconvenient for analytical treatment. We introduce /3j to decompose this matrix as 


/3j = 27r{ABJ{J + 1) - D^J\J + l)^}, 


= Ej 


= Ej 


27ri?f J(J + 1), 

7rJ(J + l) 

% ■ 


(14) 

(15) 

(16) 


aj is expressed using j3j as 


Oij = /3j+i - 13 j. 


(17) 


Using /3j, eq. ([12]) can be decomposed as 


where 


Hit) = -^Tit)Vfit), 


T{t) 


Q^oit 0 • • • 0 

0 ••• 0 

: : 0 

0 0 0 


0 1 0 


1 0 1 


V = 


0 1 0 


0 

0 

0 


1 


0 0 0 1 0 


(18) 


(19) 


( 20 ) 


where T'(t) is the conjugate matrix of T'(t), namely, T(t) •T(t) = I, where I is the identity matrix. 
Up until this point, the interaction matrix Hit) was decomposed into time-dependent diagonal 
matrices and a time-independent tridiagonal matrix. 

The matrix in ea. (|18n can be decomposed again. Using the method called spectral decomposi¬ 
tion, ea. fpn)) can be decomposed as 


V = YAY, 


(21) 
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where A is a diagonal matrix given by 

Ao 0 ••• 0 

0 Ai ••• 0 

A = 

: : 0 
0 0 0 Ajv 

Afc = 2 cos , (0 < A: < A), 

and is a unitary matrix given by 



2/00 

2/10 •• 

■ ym 

2/01 

2/11 •• 

■ ym 

yoN 

yiN •• 

■ yNN 


Vjk 




/ (j + l)(fc + l)7r \ 

V A' + l )' 


From the result above, the interaction matrix can be given as 


H{t) = -^^T{t)YAYf{t), 


( 22 ) 


(23) 


(24) 


(25) 


(26) 


D. Analytical solution with an assumption of time-order invariance 

From here, we derive an analytical solution for the time-evolution of the rotational probability 
amplitude. Generally, for the case of matrices such as eq. m, we have to consider the time-ordered 
product, which is the integral of the interaction Hamiltonian along all possible quantum paths. The 
exact solution is given by the formula as 

°° ft ftn-l ftl 

C{t) = C{to) + ^{-iT / dtn-l / dtn-2- • • / dtoH{tn-l)H{tn-2) ' ' ' H{to)C{to), (27) 

n=l ''to ''to •'to 

where C{t) is a vector representation of Cj{t). However, in our case, because many quantum states 
should be considered, infinite quantum paths must be examined. Such a mathematically exact 
treatment is impossible for the physical model presented. Here we derive a tentative expression 
for the probability amplitude by ignoring mathematical precision. We assume the time-order 
invariance of the Hamiltonian, and apply the formula, 

C{t) = C(to) + / dtn-l / dtn-2- • • / dtoH{tn-l)H(tn-2) ' ' ' H{to)C{to), (28) 

•'to •'to •'to 

= H{t)dt\ C{to). (29) 
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Using ea. (l26p . we can easily derive a tentative expression of the probability amplitude from J = r 
to J = s, given by 

OO N ^ y . '\^ { '\ ^ 

o.(*) = EEst -V / ex.p{Psit) ■ ex.p{-l3rit)dtj yskKVrk-, (30) 

n=0 fc=0 \ ' \dto ) 

OO N , „ \ 1 

“EEbra — (exp{(/3s - Pr)it})"' yskKVrk- (31) 

n=0fc=0 ^ 

An interesting characteristic of this expression is that it converges to the mathematically exact 
solution under a fully resonant condition as 


OO N 

lim Crsit) = 

n=0 fc=0 


/ iy-BftAfc \ ’ 


n! V 


VskUrk' 


(32) 


The expression in ea. (l32p is regarded as the Bessel function if the network boundaries are located 
far from the initial state. The expression derived here takes an extended form of the exact solution 
under a fully resonant condition. However, the role of 7 is changed between the fully resonant 
condition and others. In ea. (p2]) . 7 and t can be swapped, and the role of 7 is regarded as the 
velocity of wave packet propagation. However, in ea. (l30p . 7 and t cannot be swapped. 7 can be 
considered to relate directly to the probability amplitude of each state. 


E. Unified parameter 

The expression obtained in ea. (l30P cannot be used for calculating time-evolution because of 
its poor convergence property. However, it seems natural that a term in the series might have 
information relating to the amplitudes. We assume the expression below as representing the unified 
localization parameter, that is, 

7 • Bf 
\/3s-/3r 

The form of this parameter, particularly the dependence of ul on 7, / 3 r, and f3s agrees with 
our intuition for localization. This parameter is actually very useful for predicting the range of 
localization. By comparing with numerical simulations, we obtain an empirical rule for localization 
using ul. The localization is found to occur in the region 

UL > 0.5. (34) 

Outside this region, the probability amplitude decreases very rapidly. The details of validation and 
discussion will be given in the remainder of this article. 







10 


III. NUMERICAL EXPERIMENT 


Here we prove the validity of the value of obtained ul by comparison with numerical simulations. 
We mainly solve eq. m using the fourth-order Runge-Kutta method. Furthermore, we solve the 
set of eq. © and eq.dS]) when evaluating the validity of the rotating-wave approximation. We 
perform a numerical demonstration of the isotope-selective rotational excitation between 
and as an example. For the molecular constants, D^/Bm = 8.21 x 10“^ for and 

8.45 X 10-’^ for ^Sk^ScI are used [351 . The ratio Bm of ^Sr^ScI to that of ^Sr^^CI is 1.03. The 


interval of pulses is tuned to synchronize the rotational period of ^^R^^Cl. The initial rotational 
state is set to J = 5 to clearly show the localization with its upper and lower boundaries. The 
results are given as relative values of a time-averaged probability distribution for each rotational 
state after irradiation with a great many pulses. The maximum probability is always obtained for 
the initial state (J = 5), which is used as the standard for normalization. 

The number of pulses for averaging is set to by considering the pulse of 7 = 1 increments 
J by 1; the increment is proportional to 7 . Under the rotating-wave approximation, the number of 
calculation points between the two pulses is set to 10 • 7 for a resonant isotope or 20 • 7 for a non¬ 
resonant isotope. In the case without the rotating-wave approximation, the number of calculation 
points is set to 4000 • 7 or 8000 • 7 because the change in the electric field is steeper than that in 
the case with the rotating-wave approximation. The numerical errors in the sum of the rotational 
distribution are less than 10 “® in the case with the rotating-wave approximation, or less than 10 “^ 
in the case without the rotating-wave approximation. 


IV. RESULTS 

A. Dependence on 7 

Figure [3] shows the results for the time-averaged probability distribution obtained with the 
numerical calculations using eq. m with various values of 7 . The pulse interval Tp was tuned to 
be Tm of ^®R^^C1. The boundaries of the time-averaged distribution were in complete agreement 
with the line ul = 0.5, and the dependence of the range on 7 was successfully explained. With the 
rotating-wave approximation, ea. (l33p was able to completely characterize the localization of the 
distribution. 

The result of the numerical calculations without the rotating-wave approximation are shown in 
Fig. H]for comparison. The calculated time-averaged distribution begins to show small differences 
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FIG. 3. Time-averaged relative probability distribution with various 7 values calculated using eq. (ED and 
the boundary line of ul = 0.5 of (a) and (b). The interval of pulses was synchronized with 

Tyi for The color scale is the same in both (a) and (b). 


from Fig. 3 around 7 = 10, and widely disperses at 7 > 40. The unified parameter rip was still 
useful in the region 7 < 20 . 




FIG. 4. Same as Fig. [3] but calculated without the rotating-wave approximation. The scale of color is the 
same as that for Fig. [31 


As shown from the results above, up can be used to predict the upper and lower boundaries of 
the rotational distribution under the conditions of the rotating-wave approximation, regardless of 
variations in the pulse interval and centrifugal distortions. At higher 7 values, where the rotating- 
wave approximation cannot be applied, the probability distribution spreads over a wider range 
than predicted with ul in most cases. 
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B. Dependence on Tp 

Next, we evaluated the validity of ul with various values of Tp by fixing 7 . Figure [5] shows the 
results of 7 = 1 and 7 = 5 . The rotating-wave approximation was applied because the values of 7 
were sufficiently small. By extending the pulse interval slightly from the rotational period of the 
molecules, the range of the time-averaged rotational distribution became slightly wider. This is 
because slight shifts by Djj{J + 1)^ were subtly compensated by tuning the peaks of the frequency 
comb. From Fig. El^a) and (b), we can see that the best pulse interval for maximum compensation 
depends on 7 . The ul successfully represents the limit of the compensation for each 7 . 

1.0020 
1.0015 
1.0010 
^ 1.0005 

1.0000 
0.9995 
0.9990 

0 10 20 30 40 50 60 0 20 40 60 80 

J J 

FIG. 5. Time-averaged relative probability distribution with various Tp values calculated using eq. (HU and 
the boundary line of up = 0.5 of with 7 = 1(a) and 7 = 5(b). The scale of color is the same as that 

in Fig. [21 Note that both the horizontal and vertical axes are different from each other in (a) and (b). 

However, the up boundary lines deviate from the calculated distribution at Tp = 1.0005 — 1.0010 
in (a), and Tp = 1.0010 — 1.0025 in (b). The intuitive reason for this is that the local minimum of 
ul prevented wave packet propagation. In other words, for the frequency comb, the frequency gap 
between absorption and the comb becomes large in the middle region of J. 

V. DISCUSSION 

First, we made a drastic assumption by ignoring the dependence of the interaction on quantum 
paths. The analytical expression obtained could not be used for calculation of the exact time- 
evolution; however, its coefficient was observed to be a useful parameter. The reason we could 
obtain approximate information from the coefficient probably relates to the characteristics of the 
spectral decomposition. In the spectral decomposition, the Hamiltonian matrix, which depends 
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on quantum paths, is expanded into repeating diagonal interactions with a basis transformation 
at the beginning and the end. In the case without localization, the approximate time-evolution 
can be calculated using the same analytical method by ignoring the dependence on quantum paths 
[^ . In the case with localization, we probably have to obtain a first-order approximation as a 
coefficient that can be related to the approximate probability amplitude. 

There seems to be a presupposition that the unified parameter ml can be used only under 
conditions with the existence of quantum paths from the initial state to the final state. A good 
counter-example to this presupposition is presented in Fig. 5. In these calculations, there are 
isolated regions where ml > 0.5 in the upper right in both figures. If the value of ml alone 
determined the probability amplitude, this amplitude would also be observed in the isolated region. 
However, there is actually no quantum path to the isolated region from the initial state, and no 
probability is observed in the isolated region. 

A natural question to ask is why is the boundary value 0.5? In our study, unfortunately, 
this value was obtained only through empirical investigations. However, there will be a firm 
mathematical reason why that boundary value should be 0.5. 

In this study, we assume a constant transition moment between the rotational states. At M = 0, 
this assumption has little influence on the evolution of the rotational distribution. At low M, the 
unified parameter is still usable because the component of the parameter is not directly related to 
the matrix of the transition moment. Almost the same mathematical treatment can be applied 
because the transition moment matrix lea l 20 p can be always numerically diagonalized. At high 
M, the validity of the unihed parameter should be evaluated in each case because the change in 
the matrix lea l20l) can be regarded as an increase or decrease in 7 . It is likely that a more general 
theoretical treatment including the difference in M can be found. 


VI. CONCLUSION 

In this study, we derived a unified parameter to predict the range of localization of the rotational 
distribution of diatomic molecules in a train of optical pulses. The validity of the obtained param¬ 
eter was evaluated by comparison with numerical simulations. Although the theoretical treatment 
was partly tentative, the obtained parameter can almost completely explain the localization induced 
by both interval mismatch and centrifugal distortion, known as Anderson localization and Bloch 
oscillation, respectively. The odd dependence on the intensity of the pulses was also explained. 
Because the unihed parameter comprises molecular energies and simple experimental parameters, 
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the parameter is useful for the planning of experiments and development. In this study, we only 
deal with the resonant pulses in the terahertz region; however, our method can be extended to 
non-resonant pulses. The parameter derived for non-resonant pulses will predict the position of 
Jr, which was mentioned in a previous work by Flofi et al 331. 
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